#include "Quad9NodeShpFunc.h"
#include "DataStructure/Common/Real.h"

using SG::DataStructure::Common::Real;
using SG::Algebra::Matrixd;
using SG::FEMSolver::Element::Quad9NodeShpFunc;

Quad9NodeShpFunc::Quad9NodeShpFunc ()
:NSet (), DNDxiSet ()
{}

/** @brief 计算积分点处形函数及形函数对参数坐标导数
 * 
 * 缓存各个积分点处的计算结果，避免下一次积分点循环时重复计算（内存换取计算时间）
 */
void Quad9NodeShpFunc::Compute (const SG::DataStructure::FEM::IntegPointArray& points)
{
    const auto nGauss = points.size ();
    NSet.resize (nGauss);
    DNDxiSet.resize (nGauss);
    for (std::size_t iGauss = 0; iGauss < nGauss; ++iGauss)
    {
        const Real& xi = points[iGauss].m_x;
        const Real& eta = points[iGauss].m_y;

        NSet[iGauss] = Matrixd (1, 9);
        NSet[iGauss] (0, 8) = (1.0 - xi * xi) * (1.0 - eta * eta);
        NSet[iGauss] (0, 7) = (1.0 - xi) * (1.0 - eta * eta) * 0.5 - NSet[iGauss] (0, 8) * 0.5;
        NSet[iGauss] (0, 6) = (1.0 - xi * xi) * (1.0 + eta) * 0.5 - NSet[iGauss] (0, 8) * 0.5;
        NSet[iGauss] (0, 5) = (1.0 + xi) * (1.0 - eta * eta) * 0.5 - NSet[iGauss] (0, 8) * 0.5;
        NSet[iGauss] (0, 4) = (1.0 - xi * xi) * (1.0 - eta) * 0.5 - NSet[iGauss] (0, 8) * 0.5;
        NSet[iGauss] (0, 3) = (1.0 - xi) * (1.0 + eta) * 0.25 - (NSet[iGauss] (0, 6) + NSet[iGauss] (0, 7)) * 0.5 - NSet[iGauss] (0, 8) * 0.25;
        NSet[iGauss] (0, 2) = (1.0 + xi) * (1.0 + eta) * 0.25 - (NSet[iGauss] (0, 5) + NSet[iGauss] (0, 6)) * 0.5 - NSet[iGauss] (0, 8) * 0.25;
        NSet[iGauss] (0, 1) = (1.0 + xi) * (1.0 - eta) * 0.25 - (NSet[iGauss] (0, 4) + NSet[iGauss] (0, 5)) * 0.5 - NSet[iGauss] (0, 8) * 0.25;
        NSet[iGauss] (0, 0) = (1.0 - xi) * (1.0 - eta) * 0.25 - (NSet[iGauss] (0, 7) + NSet[iGauss] (0, 4)) * 0.5 - NSet[iGauss] (0, 8) * 0.25;

        DNDxiSet[iGauss] = Matrixd (2, 9);
        DNDxiSet[iGauss] (0, 8) = -2.0 * xi * (1.0 - eta * eta);
        DNDxiSet[iGauss] (0, 7) = -(1.0 - eta * eta) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.5;
        DNDxiSet[iGauss] (0, 6) = -2.0 * xi * (1.0 + eta) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.5;
        DNDxiSet[iGauss] (0, 5) =  (1.0 - eta * eta) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.5;
        DNDxiSet[iGauss] (0, 4) = -2.0 * xi * (1.0 - eta) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.5;
        DNDxiSet[iGauss] (0, 3) = -(1.0 + eta) * 0.25 - (DNDxiSet[iGauss] (0, 6) + DNDxiSet[iGauss] (0, 7)) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.25;
        DNDxiSet[iGauss] (0, 2) =  (1.0 + eta) * 0.25 - (DNDxiSet[iGauss] (0, 5) + DNDxiSet[iGauss] (0, 6)) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.25;
        DNDxiSet[iGauss] (0, 1) =  (1.0 - eta) * 0.25 - (DNDxiSet[iGauss] (0, 4) + DNDxiSet[iGauss] (0, 5)) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.25;
        DNDxiSet[iGauss] (0, 0) = -(1.0 - eta) * 0.25 - (DNDxiSet[iGauss] (0, 7) + DNDxiSet[iGauss] (0, 4)) * 0.5 - DNDxiSet[iGauss] (0, 8) * 0.25;

        DNDxiSet[iGauss] (1, 8) = -(1.0 - xi * xi) * (2.0 * eta);
        DNDxiSet[iGauss] (1, 7) = -(1.0 - xi) * (2.0 * eta) * 0.5 - DNDxiSet[iGauss] (1, 8) * 0.5;
        DNDxiSet[iGauss] (1, 6) =  (1.0 - xi * xi) * 0.5 - DNDxiSet[iGauss] (1, 8) * 0.5;
        DNDxiSet[iGauss] (1, 5) = -(1.0 + xi) * (2.0 * eta) * 0.5 -DNDxiSet[iGauss] (1, 8) * 0.5;
        DNDxiSet[iGauss] (1, 4) = -(1.0 - xi * xi) * 0.5 - DNDxiSet[iGauss] (1, 8) * 0.5;
        DNDxiSet[iGauss] (1, 3) =  (1.0 - xi) * 0.25 - (DNDxiSet[iGauss] (1, 6) + DNDxiSet[iGauss] (1, 7)) * 0.5 - DNDxiSet[iGauss] (1, 8) * 0.25;
        DNDxiSet[iGauss] (1, 2) =  (1.0 + xi) * 0.25 - (DNDxiSet[iGauss] (1, 5) + DNDxiSet[iGauss] (1, 6)) * 0.5 - DNDxiSet[iGauss] (1, 8) * 0.25;
        DNDxiSet[iGauss] (1, 1) = -(1.0 + xi) * 0.25 - (DNDxiSet[iGauss] (1, 4) + DNDxiSet[iGauss] (1, 5)) * 0.5 - DNDxiSet[iGauss] (1, 8) * 0.25;
        DNDxiSet[iGauss] (1, 0) = -(1.0 - xi) * 0.25 - (DNDxiSet[iGauss] (1, 7) + DNDxiSet[iGauss] (1, 4)) * 0.5 - DNDxiSet[iGauss] (1, 8) * 0.25;
    }
}